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ABSTRACT 

We use cosmological smoothed-particle hydrodynamics (SPH) simulations of the A 
cold dark matter (CDM) model to study the abundance of damped Lyman-a absorbers 
(DLAs) in the redshift range z — Q — 4.5. We compute the cumulative DLA abun- 
dance by using the relation between DLA cross-section and the total halo mass inferred 
from the simulations. Our approach includes standard radiative cooling and heating 
with a uniform UV background, star formation, supernova feedback, as well as a phe- 
nomenological model for feedback by galactic winds. The latter allows us to examine, 
in particular, the effect of galactic outflows on the abund ance of DLAs. We employ th e 
"conservative entropy" formulation of SPH developed bv lSpringel fc HernauistI l)2002|) . 
which mitigates against the systematic overcooling that affected earlier simulations. 
In addition, we utilise a series of simulations of varying box-size and particle number 
to isolate the impact of numerical resolution on our results. 

We show that the DLA abundance was overestimated in previous studies for three 
reasons: (1) the overcooling of gas occurring with non-conservative formulations of 
SPH, (2) a lack of numerical resolution, and (3) an inadequate treatment of feedback. 
Our new results for the total neutral hydrogen mass density, DLA abundance, and 
column density distribution function all agree reasonably well with observational esti- 
mates at redshift z = 3, indicating that DLAs arise naturally from radiatively cooled 
gas in dark matter haloes that form in a ACDM universe. Our simulations suggest a 
moderate decrease in DLA abundance by roughly a factor of two from z — 4.5 to 3, 
consistent with observations. A significant decline in abundance from z = 3 to z = 1, 
followed by weak evolution from z = 1 to z = 0, is also indicated, but our low-redshift 
results need to be interpreted with caution because they are based on coarser simula- 
tions than the ones employed at high redshift. Our highest resolution simulation also 
suggests that the halo mass-scale below which DLAs do not exist is slightly above 
10* /i~^M0 at z = 3 — 4, somewhat lower than previously estimated. 

Key words: cosmology: theory - galaxies: evolution - galaxies: formation - methods: 
numerical. 



1 INTRODUCTION 

Damped Lyman-a absorbers, historically defined as quasar 
absorption systems wit h neutral hydrogen column density 
A'hi > 2 X 10 cm~ llWolfe et al.lll986D . are one of the 
best probes of structure formation in the early universe. 
Since DLAs are dense concentrations of gas often found 
at z > 3, it is natural to suppose that they are closely 
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linked to the formation of galaxies and stars at high red- 
shift. It has become clear in recent years from the study of 
Lyman-break galaxies at z ^ 3 — 4 (e.g. lAdelberger et all 
119981: ISteidel et al.lll99l IShaoley et alJl200ll) that the as- 
sembly of galaxies is actively going on at z ~ 3, consis- 
tent with hierarchical structure formation in a cold dark 
matter universe (e.g. ^ o & Fukugita 1996; Baugh ct 
[1998': 'Jin g fc SutdlT 998: Katz. Hornauist & Weinberg iggj" 
i Kauffmann et al.lll99a:.Mo. Mao. &: Whit&.199a: .Nagaminj 
120021: IWeinberg. Hernquist fc Katdl2002^ . 

A picture of the history of cosmic star formation 
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emerging from both theory and observation is that it 
rises from the present towards high redshift, even beyond 



z = 


3 (e.g. 


iPascarelle. Lanzetta. & Fernandez-Sotd Il998l 


Blain et al. 


'1999"; 'Nagamine et al.' 


'2001a; Lanzetta ct al.' 


200a 
2OO2I 


ISoringel & Hernauist ^2003h.; liornquist & Springel 
. The conversion of gas into stars is, therefore, taking 



place at a significant rate at z > 3. If DLAs dominate the 
neutral hydrogen gas content of the Universe at z ~ 3, they 
are thus serving as an important reservoir of neutral gas for 
star formation. Determining the physical nature and number 
density of DLAs may hence be one of the most important 
keys for further constraining the cosmic star formation his- 
tory and theories of galaxy formation. 

Although the current sample of observed DLAs at 2 > 
1.5 is not yet as large as that of Lyman-break galaxies (where 
« 1000 are known), the total number of DLAs that have 
been discovered is now approaching « 100, and the number 
density per unit redshift of hig h column density systems 
appe ars to peak at around z ~ 3 iStorrie-Lombardi fc Wolfa 
1200(3) . At lower redshift, the situation is rather different. The 
identification of DLAs a.t z < 1.5 has been difficult because 
they are rare and the need for ultra-violet (UV) spectroscopy 
to detect Lya absorption at these low redshifts. The number 
of quasars studied in UV has been smal l until recently. To 
overcome this difficultv. lRao fc Turnshe^ i2000h searched for 
DLAs in 87 Mg 11 absorbers, and uncovere d 12 new systems. 
There are 23 DLAs &t z < 1.65 listed bv lRao fc Turnshekl 
J2000l) . 

Despite the likely importance of DLAs and the accumu- 
lation of observational data on them, their true nature re- 
mains controversial. Historically, it has often been suggested 
that high-redshift DLAs are large, rapidly rotating discs, be- 
cause DLAs have properties similar to local galactic discs, 
such as large neutral hydrogen column densities together 
with low degrees of ionisation and small velocity dis per- 
sions ( Wolfe et al. 1986). More recently, ^rochaska fc Wol^ 
(Il997i Il998i) argued that the observed distribution of veloc- 
ity widths and the asymmetric absorption profiles of low- 
ionisation ionic species can be best described by massive, 
rapidly rotating cold discs. 

O n the other hand, iHaehnelt. Steinmetz. fc RauchI 
lll998h examined a small number of dark matter haloes in 
a high-resolution (sub-kpc) SPH simulation, and showed 
that such observational signatures can also be explained 
by a mixture of rotation, random motions, infall, and 
mergers of protogal actic clumps. There are some obser- 
vational indications dLe Brun et al.lll997l : iRao fc Turnshe^ 
ll998l:lKulkarni et alj200ll200 jl from direct imaging studies 
that luminous disc galaxies may not represent the dominant 
population of DLA galaxies (i.e. galaxies that host DLAs). 
Although the possibility of artifacts due to point-spread- 
function effects cannot be fully excluded, these observations 
suggest that some DLA galaxies could be compact, clumpy 
objects, or low surface brightness galaxies, rather than large, 
well- formed protogalactic disks or spheroids. 

Robust numerical estimates of DLA properties have 
been hampered by the significant requirements on numerical 
resolution needed to capture t he full pop ulation of DLAs. 
Earlie r studies by iKatz et alj (Il996'l and iHernauist et alJ 
(11996?) showed that the observed Hi column density distribu- 
tion can be reproduced within a factor of a few in hydrody- 
namic simulations based on a CDM model over a wide range 



of column densities 10^*cm~^ < Nhi < lO'^'^cm"^. Their re- 
sults demonstrated that the Ly-a forest develops naturally 
in the hierarchical clustering scenario of CDM universes, 
and that DLAs and Lyman-limit systems (A^hi > lO^^cm"^) 
arise in these models from radiatively cooled gas inside dark 
matter haloes that host forming galaxies at high redshift. 
However, their calculations were based on simulations of a 
critical-density universe with Q.^ = 1, and could not resolve 
haloes with masses 

Subsequently, iGardner et"ai] (Il997allhl |200ll) extended 
the ea rlier results of lKatr^^Lr(u'996ft and Memouist et alJ 
(Il996l) by developing a method to correct for the resolution 
limitations of the simulations. They measured the relation 
between absorption cross-section and halo circular veloc- 
ity from hydrodynamic simulations, and then convolved it 
with t he analytic halo mass fu nction ("e.g. lPress fc Schechteil 
I1974I: ISheth fc TormenI Il999^ to compute the cumulative 
abundance of DLAs. Using this correction method, they 
were able to reproduce the observed abundance of DLAs 
if they required that haloes with circular velocity Vc < 
60kms"^ (which corresponds to Af Ri 2 x lO"' fe'^M©) did 
not harbour DLAs. However, their simulations could not re- 
solve haloes with masses below 10^" /i^^Mq, although such 
haloes may still host a significant number of DLAs. If the ab- 
sorption cross-section of haloes with M < 10^° h~^MQ does 
not follow the same relation between the cross-section and 
the halo mass as determined from higher mass haloes, the 
DLA abundance could be either over- or underestimated. 
Simulations of higher resolution are hence clearly needed to 
make more robust predictions for the DLA abundance at 
redshift 2 = 3. 

Recently, ISpringel fc Hernauist! (|20o2) developed a 
novel formulation of SPH that is based on integrating t he en- 
tropy as an independe nt thermodynamic variable (e.g. lLuc'vl 
119771: iHernauistilTgg^) . and which takes variations of the 
SPH smoothing lengths self-consistently into account. They 
showed that this new version maintains contact discontinu- 
ities (as they arise at the interface between cold and hot gas 
in haloes) much better than previous treatments of SPH. 
Consequently, their formulation does not suffer from the se- 
vere overcooling that was typically seen in previous SPH 
simulations.^ We hence use this new methodology for our 
studies of DLAs. 

Our simulations also include a novel method for 
tre ating star format io n and feedback, as proposed by 
iSpringel fc Hernauist! (!2003a!) . It is based on a sub- 
resolution multi-phase description of the dense, star-forming 
interstellar medium (ISM) and a phenomenological model 
for strong feedback by galactic winds. The inclusion of 
winds was motivated by t he possibility that o utfiows from 
galaxies at high redshift jPettini et alj !2002!) play a role 
in distributing meta ls into the intergalactic medium (e.g. 
!Aguirre et alj2001a!lbl) . and they may also alter the distribu- 
tion of neutral gas around galaxies (jAdolbcrgcr ct al. 2003), 
although this process remains uncertain (e.g. ICroft et alJ 



^ The likelihood that earlier SPH studies were affected by over- 
cooling due to numerical effects is supported by comparisons be- 
tween our new formulation and simulations using an adaptive 
mesh refinement (AMR) algorithm (M. Norman, private commu- 
nication). This finding will be presented in due course. 
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Run 


BoxsizG 


N 




^gas 








R3 


3.375 


2 X 1443 


9.29 X 10^ 


1.43 X 10^ 


0.94 


4.00 


strong 


R4 


3.375 


2 X 216^ 


2.75 X 10^ 


4.24 X 10** 


0.63 


4.00 


strong 


03 


10.00 


2 X 1443 


2.42 X 10^ 


3.72 X 10*5 


2.78 


2.75 


none 


P3 


10.00 


2 X 144^ 


2.42 X 10'^ 


3.72 X 10^ 


2.78 


2.75 


weak 


Q3 


10.00 


2 X 1443 


2.42 X 10'^ 


3.72 X 10^ 


2.78 


2.75 


strong 


Q4 


10.00 


2 X 216^ 


7.16 X 10** 


1.10 X 10** 


1.85 


2.75 


strong 


Q5 


10.00 


2 X 3243 


2.12 X lO'' 


3.26 X 10^ 


1.23 


2.75 


strong 


D4 


33.75 


2 X 216^ 


2.75 X 10** 


4.24 X 10^ 


6.25 


1.00 


strong 


D5 


33.75 


2 X 3243 


8.15 X 10^ 


1.26 X 10^ 


4.17 


1.00 


strong 


G4 


100.0 


2 X 216-^ 


7.16 X 10^ 


1.10 X lO'' 


12.0 


0.00 


strong 


G5 


100.0 


2 X 3243 


2.12 X 10^ 


3.26 X 10* 


8.00 


0.00 


strong 



Table 1. Simulations employed in this study. The box-size is given in units of Mpc, A'^p is the particle number of dark matter and 
gas (hence X 2), moM and rrigas are the masses of dark matter and gas particles in units of /i~^Mq, respectively, e is the comoving 
gravitational softening length in units of kpc, and Zqikj is the ending redshift of the simulation. The value of e is a measure of spatial 
resolution. From the top to the bottom row, we refer to R3 & R4 collectively as 'R-series', the next 5 runs (03 to Q5) are called 'Q-series', 
D4 &: D5 are called 'D-series ', and G4 & G5 are called 'G-series'. The 'strong-wind' simulations form a subset of the runs analysed by 
ISoringel fc Hernauisj 



l2002l: iKollmeier et aljliool . Together with the increase in 
numerical resolution provided by our simulations, it is of in- 
terest to see how refinements in physical modelling modify 
the predictions of DLA properties in a CDM universe. 

In this paper, we focus on the abundance of DLAs in the 
redshift range z = — 4.5. The present work extends and 
co mplements earlier num erical work by ,Katz ct al. ( 1 99^ 
and lGardner et all i200l[) . Physical properties of DLAs such 
as their star formation rates, metallicities, and their relation 
to galax;ies will be presented elsewhere. 

The paper is organised as follows. In Section |5| we 
briefly describe the numerical parameters of our simulation 
set. We then present the evolution of the total neutral hy- 
drogen mass density in the simulations in Section |3 In Sec- 
tion |11 we describe how we compute the Hi column density 
and DLA cross-section as a function of total halo mass. In 
Section]^ we determine the cumulative abundance of DLAs, 
and discuss the evolution of DLA abundance from z = 4.5 
to 2 = 0. The Hi column density distribution function is 
presented in Section |S| Finally, we summarise and discuss 
the implication of our work in Section Q 



2 SIMULATIONS 

We analyse a large set of cosmological SPH simulations that 
differ in box size, mass resolution and feedback strength, 
as summarised in Table In particular, we consider box 
sizes ranging from 3.375 to 100/i~^Mpc on a side, with 
particle numbers between 2 x 144'' and 2 x 324^, allow- 
ing us to probe gaseous mass resolutions in the range 
4.2 X lO"* to 1.1 X 10^ h~^MQ. These simulations are partly 
taken from a study of the cos mic star formation history by 
ISprineel fc Hernauisd J2003bl) . supplemented by additional 
runs with weaker or no galactic winds. The joint analysis of 
this series of simulations allows us to significantly broaden 
the range of spatial and mass-scales that we can probe com- 
pared to what is presently attainable within a single simu- 
lation. 

There are three main new features to our simulations. 



First, we use a new "conservative entropy" formulation of 
SPH ( Springcl & Hcrnquist 2002) which explicitly conserves 
entropy (in regions without shocks), as well as momentum 
and energy, even when one allows for fully adaptive smooth- 
ing lengths. This formulation moderates the overcooling 
problem present in earlier formulations of SPH (sec also 
O^'oshida ct al. 2002; Pcarcc ct al. 1999; Croft ct al. 2001). 

Second, highly over-dense gas particles are treated with 
an effective sub-resolution model for the ISM, as described 
by SorinEcl & Hcrnauist (200^). In this model, the dense 
ISM is pictured to be a two-phase fluid consisting of cold 
clouds in pressure equilibrium with a hot ambient phase. 
Each gas particle represents a statistical mixture of these 
phases. Cold clouds grow by radiative cooling out of the hot 
medium, and this material forms the reservoir of baryons 
available for star formation. Once star formation occurs, the 
resulting supernova explosions deposit energy into the hot 
gas, heating it, and also evaporate cold clouds, transferring 
cold gas back into the ambient phase. This establishes a tight 
self-regulation mechanism for star formation in the ISM. 

Third, we implemented a phenomenological model for 
galactic winds in order to study the effect of outflows on 
DLAs, galaxies, and the intergalactic medium (IGM). In 
this model, gas particles are stochastically driven out of the 
dense star-forming medium by assigning an extra momen- 
tum in random directions, with a rate and magnitude chosen 
to reproduce mass-loads and wind speeds similar to those 
observed. See Spri ngel fc Hernguist (■2003al for a detailed 
discussion of this method. 

Most of our simulations employ a "strong" wind of 
speed 484kms"\ but for the 10/i"^Mpc box (runs 03, P3, 
Q3, Q4, Q5; collectively called 'Q-series') we also varied the 
wind strength. Therefore, this Q-series can be used to study 
both the effect of numerical resolution and the consequences 
of feedback from galactic winds. The runs in the other sim- 
ulation series then allow the extension of the strong wind 
results to smaller scales ('R-Series'), or to larger box-sizes 
and hence lower redshift ('D-' and 'G-Series'). Our naming 
convention is such that runs of the same model (box-size and 
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included physics) are designated with the same letter, with 
an additional number specifying the particle resolution. 

Our calculations include a uniform UV back g round 
radiation field of a modified iHaardt fc Madaul lll996ll 
spec trum, where rei onisation takes place at 2; ~ 6 
see iDave et alJll999^ . as suggested by observations (e.g. 
iBecker e^Lll200l[) and radiative transfer calculations of the 
imp act of the stellar sour ces in our simulations on the IGM 
fe.g. lSokasian et al. '200?). The ra diative co oling and heating 
rate is computed as described bv lKatz et al.. (■199&) assum- 
ing that the gas is optically thin and in ionization equilib- 
rium. The abundance of different ionic species, including H°, 
He°, H"*", He"*", and He"*"^, is computed by solving the net- 
work of equilibrium equations self-consistently with a speci- 
fied UV background radiation. The results presented in this 
paper should not be affected by the assumption of ioniza- 
tion equilibrium as we are dealing with high density regions 
where this assumption is satisfied well. The adopted cos- 
mological parameters of all runs are (f2m, f2A, fib, erg, /i) = 
(0.3,0.7,0.04,0.9,0.7). The simulations were performed on 
the Athlon-MP cluster at the Center for Parallel Astrophys- 
ical Computing (CPAC) at the Harvard- Smithsonian Cen- 
ter for Astroph ysics, using a modified version of t he parallel 
GADGET code iSpringel. Yoshida fc Whitdl200ll) . 



3 NEUTRAL HYDROGEN MASS DENSITY 

In our simulation methodology, gas is subject to a thermal 
instability yielding a multiphase medium once the physical 
gas density lies above a threshold density pth, which marks 
the onset of cold cloud formation. If the physical density is 
lower than this threshold, a particle represents ordinary gas 
in a single phase. In this latter case, the neutral hydrogen 
mass of the particle can be computed as follows: 

mm = A^h ■ Xh ■ mgas (p < pth), (1) 

where Xh — 0.76 is the primordial mass fraction of hydro- 
gen, and TVh is the number density of neutral hydrogen atoms 
in units of the total number density of hydrogen nuclei. The 
quantity A^h is computed by solving the ionisation balance 
as a function of density and current UV background flux. 

If the gas density is higher than the threshold density, 
we identify the mass of neutral gas with the mass of cold 
clouds contained in the multiphase medium. The mass of 
neutral hydrogen of such a multiphase particle is then given 

by 

mm = a; • Xh • mgas (p > Pth), (2) 

where x is the mass fraction of cold clouds. In the multiphase 
model of Springel & Hernquist (2002b), x can be computed 
as 

--P=/P = l + ^-/^ + ^> (3) 
where the quantities y and pth are defined by 

^" p[f3usN - {1 - PW] ^ ' 
and 

_ Xth (3USN - (1 - I3)uc , , 

" (1 - xth)^ f2A(usNMo) ■ ^ ^ 



Here, Anct(p, it) is the usual cooling function, and we have 
defined A(p,u) = Anet(p, m)/p^. The parameter /3 gives the 
mass fraction of short-lived stars that explode as supernovae, 
usN describes the energy released by the supernovae, Uc is 
the assumed temperature of the cold clouds, Uh the temper- 
ature of the hot medium, Ao the cloud evaporation parame- 
ter, and tl gives the star formation time-scale. The quantity 
Xth = {uh — U4,)/{uh — Uc) ~ 1 — AoUi/us-N is the mass 
fraction in cold clouds at the threshold density, where U4 is 
the specific energy corresponding to T = IC* K. We refer to 
Springel & Hernquist (2002b) for a more detailed explana- 
tion of these parameters, and a derivation of equations 
to 

In Figure we show the total neutral hydrogen mass 
density as a function of redshift. The values plotted are given 
in terms of flui x 10^. At redshifts above six, essentially all 
the hydrogen in the simulation box is still neutral {Qm — 
Xu^h = 0.0304). Once the ionising background sets in at 
= 6, the neutral hydrogen starts to become ionised and 
the neutral fraction decreases rapidly by a few orders of 
magnitude; i.e. reionisation takes place. Note that with the 
exception of the R-Series, our earliest simulation output at 
2 < 6 corresponds to z — 5.5. This is why most of the 
models in Figure seem to show a rapid rise of ^hi already 
at 2 = 5.5, which simply arises by drawing a line to the data 
point at z — 6 which lies a few orders of magnitude higher 
in this linear plot. 

We also include observatio nal data points, for com- 
pariso n. The data points from IStorrie-Lombardi fc Wol^ 
(|200Cf, open squares) account only for DLAs, but those 
of Pcroux ot al. (2001, filled triangles) include a correc- 
tion for the neutral gas that is in sub-PL As (10^^ < 
A^Hi < lO^^cm-^). Data points froi nlRao fc Turnshekl (l2000l. 
solid triangles) at low-redshift and lZwaan et al .1 lll99l open 
cross) at z = are also shown. 

In the left panel of Figure we compare results oidy 
for the Q-series (10 /i~'^Mpc box), allowing us to assess con- 
vergence as a function of mass resolution and to investigate 
the dependence of the results on the strength of feedback 
from winds. 

The comparison between Q3, Q4, and Q5 (strong wind) 
shows that there is quite good agreement between runs with 
different numerical resolution. In fact, the results for Q3, Q4, 
and Q5 are essentially identical at z ~ 4.5, with Q5 being 
slightly higher for z > 4.5 than its lower resolution counter- 
parts, while for z < 4.5 the opposite trend is observed. This 
mild effect can be understood as follows: in a higher resolu- 
tion run like Q5, many more small dark matter haloes can 
be resolved than in a lower resolution run (like Q3) , partic- 
ularly at early times, where gas can cool very efficiently. As 
a result, the higher resolution simulation develops a larger 
fraction of cold and hence neutral gas at high redshift. How- 
ever, an increased fraction of cold gas will also trigger more 
intense star formation that both consumes neutral gas and 
leads to gas ejection by winds from low mass haloes. Subse- 
quently, the neutral fraction can then fall slightly below the 
lower resolution runs. 

Comparing the results for the models 03 (no wind), 
P3 (weak wind), and Q3 (strong wind) shows the impact of 
feedback by galactic winds. As the wind strength increases, 
the neutral density flm decreases. More neutral gas is then 
ejected out of the dense ISM into the intergalactic medium. 
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Q-Series 



a 



n \ \ I \ r 

"03 (no wind) 

-P3 (weak wind) 

-Q3 (strong wind) 

-Q4 (strong wind) 

-Q5 (strong wind) 



- * 





2 4 6 

redshift 

Figure 1. Evolution of the total neutral hydrogen mass density in each simulation box as a function of redshift. The plot ted values are 
Qht X 10^. We also show observational data points from Istorrie- Lombardi fc Wolfd <200l open squares; on l y for DLAs"). |Peroux et alj 
I2OOII, filled t riangl es; including the correction for the neutral gas not included in DLAsV lRao fc Turnshekl i200(l open triangles), and 
IZwaan et al.l il997t open cross at 2; = 0). Left panel: A comparison of the Q-series (runs in lO/i^'^Mpc boxes) is shown. The decrease in 
Ojji from 03 (no wind) to P3 (weak wind), and then to Q3 (strong wind), shows the effect of feedback by galactic winds. The comparison 
between Q3, Q4, and Q5 shows the level of convergence achieved for runs with different resolution. For P3 (top short-dashed line), we 
separately show Hhi in regions of overdensity 1 + 5 > 10'^ and 10* (middle and bottom short-dashed line, respectively). Right panel: 
Results for the R-, D-, and G-series are shown. Q5 is also included to ease comparison with the left panel. Results for R3 and G4 are 
omitted for clarity (see text). 



where it becomes highly ionised by UV background radia- 
tion. Interestingly, '03' (no wind run) exceeds all observed 
data points, so a feedback eflect such as galactic winds ap- 
pears necessary to make the f^ni measurements of the sim- 
ulations consistent with observations. The results for our 
'strong-wind' runs (Q3, Q4, Q5) underpredict the observa- 
tional estimates at z = 3 slightly, but there is still marginal 
agreement within 1 cr, which is encouraging. However, the 
best value for the galactic wind strength parameter for our 
simulation seems to lie somewhere between that of P3 (weak 
wind) and the Q-runs (strong wind). 

For the 'P3' run, we also show separate measurements 
of fiHi restricted to regions of overdensity 1 -\- 5 > 10^ and 
10*, respectively (red short-dashed lines). The fact that the 
lines for 1 + <5 > 10* and 10^ have converged by z ~ 3 shows 
that most of the neutral hydrogen mass in the universe is 
already in a highly concentrated form by this epoch. 

In the right panel of Figure we show our results for 
simulations of the R-, D-, and G-series, together with Q5 
for reference to the left panel. The results for D4 and D5 
are consistent with one another at 2 = 3. 'R3' is not shown 
because it is almost identical to 'R4', and 'G4' is omitted 
because it underpredicts Qhi significantly due to lack of res- 
olution at z > 3. By comparing to the simulations of the Q- 
and D-series, we see that the resolution of the G-series is not 
sufficient to correctly describe the neutral fraction at z = 3. 
This is because even the 2 x 324^ run G5 misses the neutral 
gas content in large numbers of small dark matter haloes 



that are present in the higher resolution runs at z = 3, such 
as those of the Q-series. Therefore, we consider Q5 to be the 
most reliable run at 2 = 3 among our simulation set. We 
also see that Q,m of 'R4' is lower than that of 'Q5', despite 
the fact that the R-series has higher mass resolution than 
the Q-series. This is likely due to the rather small box-size 
of the R-series compared to the Q-series, which leads to an 
insufficient sampling of rare, massive objects, and compro- 
mises the use of R4 as a truly representative sample of the 
universe. 

The effect of the multiphase model adopted in the cur- 
rent simulations can be assessed by setting the value of cold 
gas mass fraction to a:: = 1 for the multiphase gas parti- 
cles [see Equation ||3J]. We find that the value of Q,m be- 
comes larger by about 15% in such a case. This suggests that 
previous formulations of hydrodynamic simulations with- 
out a consideration for the multiphase nature of the gas 
would have overestimated the cold gas fraction by a similar 
amount . 



4 Hi column DENSITY & DLA 
CROSS-SECTION 

We now describe how we compute the Hi column density 
A^'hi and the DLA cross-section ctdla for each dark matter 
halo. First, we identify dark matter haloes by applying a 
conventional friends-of-friends algorithm to the dark matter 
particles in each simulation. We set the minimum number of 
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dark matter particles for a halo to 32; i.e. haloes with fewer 
particles are not included in the group catalogue. We have 
confirmed that the dark matter halo mas s functions agree 
well w ith the analytic mass function of 'Sheth fc To rmenI 
1^9^. After dark matter haloes are identified, we associate 
each gas and star particle with their nearest dark matter par- 
ticle, including them in the particle list of the corresponding 
haloes, when appropriate. 

Then, for each halo, a uniform grid covering the entire 
halo and whose grid-size is equal to the gravitational soft- 
ening length, is placed at the center-of-mass of the halo. We 
then project the neutral gas in the halo onto a plane, and 
obtain the column density of each grid-cell in this plane. 
The neutral mass of each gas particle is smoothed over a 
spherical region of grid-cells, weighted by the SPIf kernel. 
To check the robustness of the result, we also tried a cloud- 
in-cell assignment scheme where the neutral mass of each gas 
particle is uniformly distributed over a cubic region of size 
^ = {^)^^^s and i(^)^/^s centered on the particle. Here 
s is the SPH smoothing length. Differences in the smooth- 
ing method can lead to slight differences in the Hi column 
density distribution, as we will discuss later in Figure |7| In 
the following, we adopt the SPH smoothing method for our 
primary results unless explicitly stated otherwise. 

Once the comoving neutral mass density pi,m in each 
grid-cell of volume is known, it is straightforward to 
project the density distribution along the direction perpen- 
dicular to the plane to obtain the column density as 

iVm = pi,m ■ e/nip ■ (1 + z)^ , (6) 

i 

where t is the comoving gravitational softening length, rUp 
is the proton mass, and z is the redshift. 

Note that in the present study, we do not apply a self- 
shielding correction when computing the neutral hydrogen 
fraction. As JCatz ct al. ( 1996) have shown, damped systems 
with column densities above A^hi — lO^'^cm"^ are essen- 
tially fully neutral and are not affected by self-shielding. 
The correction is expected to be large for systems with 
10" < Ahi < lO^^cm"^, however. A full 3-dimensional 
treatment of self-shielding is beyond the scope of the present 
study, but it is clearly an interesting and important issue 
in its own right. The results presented in this paper for 
very high column density systems should however be robust 
against self-shielding corrections. 

Once the column density of each cell in the projected 
plane is obtained, we estimate the comoving DLA cross- 
section of each halo by simply counting the number of grid- 
cells that exceed Ahi = 2 x 10^°cm~'^ and multiplying this 
number by the comoving unit area of the grid-cells. 

4.1 DLA cross-section at redshift 3 

In Figure 121 we show the comoving DLA cross-section ctdla 
as a function of total halo mass Mtot at redshift 2 = 3. All 
panels are for runs that include 'strong winds', except where 
explicitly labeled otherwise. The data points are binned in 
terms of log Mtot , and the median value in each bin is shown 
by the open triangles. The quartiles in each mass bin are 
shown as error bars. For plotting purposes only, we assign an 
arbitrary value of log ctdla = 0.3 to all haloes with no DLAs, 
and they are shown by the crosses at the bottom of each 



Redshift 3 



Run 


slope a 


P 


03 


0.72 


3.94 


P3 


0.79 


3.99 


Q3 


0.84 


3.98 


Q4 


0.93 


4.03 


Q5 


1.02 


4.18 


D4 


0.68 


3.93 


D5 


0.81 


3.96 


G4 


0.64 


3.88 


G5 


0.62 


3.89 



Table 2. The parameters obtained by least-square fitting the 
power-law of equation Q to the median points shown in Figurel^ 



panel. We included these 'no-DLA haloes' in computing the 
median cross-section, therefore the error bars in the lowest 
mass-bins sometimes extend to the bottom of the figure. 

We then fit the median points to a power law, ctdla oc 
i\ftot", assuming a functional form of 

log CTDLA = a(log A/tot - 12) -I- /3, (7) 

and determine the values of the slope 'a' and the normalisa- 
tion '/3' by least-squares fitting. The value of /3 hence gives 
the value of log ctdla at Mtot = 10^^ /i~^M0. We chose this 
reference mass-scale because it is well covered by most of 

the simulations used in this paper^ 

Unlike the analysis of Gard ner et al.l I^QT^), we do not 
invoke a limiting halo mass below which a dark matter halo 
does not harbour a DLA. As can be seen in all the panels 
of Figure |21 such a clear cutofi' does not really exist, and 
DLAs continue to be found in haloes with masses down to 
Mtot ^ 10*-^ /i-^Mq in 'Q5'. We will come back to this point 
later. 

We summarise the results of our power-law fitting in 
Table |5| It is satisfying that the values of the normalisa- 
tion '/3' agree very well among different runs. This demon- 
strates that our results for runs with widely varying reso- 
lution are numerically well-converged at the mass-scale of 
Mtot = 10^^ /i~^M0. It is seen that the slope 'a' becomes 
steeper as the galactic wind strength increases from 03 to 
P3, and further to Q3. This is because gas in low-mass haloes 
is lost at a higher rate in runs with stronger winds, making 
the DLA cross-sections decrease for small haloes. Another 
trend seen in Table |5| is that the slope becomes somewhat 
steeper as the resolution of the simulation increases from Q3 
to Q4, and then to Q5. This can be explained by the fact 
that higher resolution simulations can resolve star formation 
in small haloes, leading to ejection of gas out of them, lower- 
ing their content of neutral gas. On the other hand, a lower 
resolution simulation misses this star formation, resulting in 
an overestimate of the baryon and neutral fraction in the 
first generation of haloes that is 'seen' in the simulation. 

iGardner et alJ i200 reported a slightly shallower slope 
even compared to '03' (no- wind run): ctdla oc v].'^"^ oc M"'^^ 
(see Table 2 of their paper). Here, Vc is the circular velocity 
of a halo, related to the halo mass by oc M"-'^. A shal- 
lower slope in general implies a higher abundance of DLAs. 
A number of effects are responsible for this difference: (1) 
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Figure 2. Comoving DLA cross-section (Tj3la. a-s a function of total halo mass at z = 3 for the Q-, D-, and G-series. All panels are for 
runs with 'strong winds', except where labelled otherwise. We show median values for logarithmic mass bins as open triangles, and the 
quartiles in each bin are indicated in the form of error bars. The horizontal concentration of crosses at log cr jjla = 0-3 indicates haloes 
without DLAs. The solid lines are power-law fits to the median points as described in the text. 
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Their resolution was slightly lower (A'^ — 128^) than that 
of our 03/P3/Q3-runs. (2) The overcooling problem in pre- 
vious formulations of SPH may have caused the slope to 
be shallower owing to an overestimate of the neutral gas 
fraction, particularly in small haloes. With our new 'conser- 
vative entropy' formulation of SPH, the cold gas fraction in 
haloes is expected to be lower, although the magnitude of 
this effect as a function of halo mass is not fully clear. (3) 
Their treatment of feedback is known to be inefficient, be- 
cause thermal energy injected into the gas is radiated away 
very rapidly. However, given that our '03'-run overpredicts 
f^Hi (Figure0l at z = 3, some form of strong feedback seems 
necessary to provide agreement with the observations. Not- 
ing that the slope of the power-law fit steepens as the wind 
strength an d resolution increase, we hence conclude that 
the slope of iGardner et all |20^) was probably too shal- 
low. This conclusion will be strengthened when we discuss 
the abundance of DLAs in Sectionj^and the column-density 
distribution function in Section |H] 



Redshift 4.5 



Run 


slope a 


/3 


R3 


0.79 


4.29 


R4 


0.86 


4.46 


03 


0.64 


4.21 


P3 


0.67 


4.24 


Q3 


0.71 


4.31 


Q4 


0.79 


4.44 


Q5 


0.87 


4.59 


D4 


0.68 


4.34 


D5 


0.74 


4.42 


G4 


0.56 


4.31 


G5 


0.65 


4.30 



Table 3. The parameters obtained by least-squares fitting the 
power-law of equation Q to the median points shown in FigurelSl 
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z=4.5 
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Figure 3. Same as Figure|5] but for z = 4.5. Here, instead of G4 and G5, we sliow results for R3 and R4. Tlie solid slanted lines are the 
best-fits to the median points (indicated by the open triangles) with a power-law f eg nation 171. and the short-dashed lines are the fits 
we obtained at z = 3, for comparison. The horizontal concentration of crosses at logo-jjLA = 0.3 indicates haloes that do not harbour a 
DLA. Note that the range along the x-axis is different from that of Figure 13 



4.2 DLA cross-section at redshift 4.5 

In Figure |3| we show the DLA cross-section a,t z = 4.5 as a 
function of total halo mass. As before, the solid lines show 
power-law fits that were obtained as described in the previ- 
ous subsection, while the short-dashed lines are the fits at 
redshift z — 3, for comparison. Note that we do not plot 
results for the G-series but instead show the R-series, which 
has much higher resolution, but was evolved only to z — 4 
due to its small box-size. The results of the power-law fit- 
ting are summarised in Table|3 Similar to z = 3, the values 
of the normalisation '/3' agree very well with each other be- 
tween runs of differing resolution. The generic trends in the 
slope as a function of wind strength and resolution are also 
similar to what we saw for z = 3. 

At the low-mass end of runs R3 and R4, we see that 
the DLA cross-section is dropping prominently at Mtot ~ 
10*'^ h~^MQ, strongly departing from the fitted power-law. 
To avoid being affected by this turn-down, we do not use 
the median points below logcroLA = 1.0 for our power- 
law fitting. A similar feature is seen in Q4 and Q5 at the 
same mass-scale. Note that the resolution of the R-series is 
sufficient to resolve haloes with Mtot = 10* h~^MQ quite 



well, so this downturn in DLA cross-section is unlikely to 
be a resolution artifact. Instead, it is probably due to the 
physical effect that the gas in these haloes is easily photoe- 
vaporated by the ionising background, and/or ejected from 
haloes due to supernovae feedback. At z = 3, the mass-scale 
of 10* — 10*'^ h~^MQ corresponds to a circular velocity of 
10 — 15 kms~^, or a virial temperature slightly below 10* K. 
Note that this is a smaller mass-sc ale th an was suggested by 
iQuinn. Katz. fc Efstathioul Jl99dl and iThoul fc Weinberd 
(T99^, who argued that haloes with circular velocities less 
than 40kms~^ are unlikely to harbour DLAs. We will dis- 
cuss this point further in Section [7| 



4.3 DLA cross-section at lower redshift 

In Figure 2] we show DLA cross-sections as a function of 
total halo mass for z = 1 and z — 0, and the parameters of 
the fitted power-laws are summarised in Table 2] A similar 
trend in the slope as a function of resolution exists at z = 1 
as we saw at z = 3. It is clear that the slope cannot be 
determined reliably for G4 and G5 at z = (and possibly 
at z = 1) due to limited resolution, as is evident from the 
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Figure 4. Same as Figure |3] but for z = 1 and 0. The short- 
dashed Unes are the fits we obtained at z = 3, for comparison. 



Lower Redshifts 



Run 


slope a 




z=l 


D4 


0.74 


2.82 


D5 


0.82 


2.75 


G4 


0.61 


2.93 


G5 


0.67 


2.89 


z=0 


G4 


0.54 


2.73 


G5 


0.56 


2.62 



Table 4. The parameters obtained by least-square fitting the 
power-law of equation to the median points shown in Figurel4l 



'stripes' seen at low cross-sections in the bottom two panels 
of Figure 3] 



5 CUMULATIVE ABUNDANCE OF DLAS 

Dark matter haloes with masses below the resolution limit 
of a simulation cannot be resolved. This is a serious prob- 
lem when one tries to compute the number density of DLAs 
based on a cosmological simulation that does not resolve all 



small mass haloes that may host a DLA. Note in particu- 
lar that the number density of dark matter haloes is known 
to increase strongly towards lower masses. Even a small in- 
completeness at low masses will hence prevent a reliable es- 
timate of the DLA abundance if only a simple number count 
of DLAs found in a cosmological s imulation is used. 

T o overcome this limitation, iGardner et al ] lll997allbl 
I2OOII) convolved a theoretical fit to the dark matter halo 
mass function with the measured relationship between DLA 
cross-section and halo mass. In this way, they were able to 
correct for incompleteness in the resolved halo abundance of 
the simulations. The cumulative abundance (or equivalently 
the rate of incidence) of DLAs per unit redshift as a function 
of halo mass in this approach can be expressed as 

— ^(> M, z) = -f nd,n(M', z) auLA{M', z) dM' , (8) 
az az J 

where nd-m{M,z) is t he dark matter halo m ass function 
(for which we use the ISheth fc TormenI 1^9^ parameter- 
isation), and Ar/Az — c/H{z) with H{z) = HoE{z) = 
Ho ^nm(l + ^)3-f for a fiat universe. In order to carry 
out this integral, the power-law fits obtained in Section 2] 
can be used to represent ctdla(M, 2) which give the mean 
relation between the halo mass and the DLA cross-section. 
Note that the dependence on the Hubble constant disap- 
pears on the right-hand-side of equation (8) because dr/d2: 
scales as , while ndmdA/ depends on , and ctdla scales 
as in the simulation. 

5.1 DLA abundance at redshift 3 

In Figure |3 we show the cumulative abundance of DLAs 
per unit redshift at 2: = 3 as a function of total halo mass. 
The horizontal shaded region in the left pan e l ind icates 
the observed DLA abundance of Pcroux et al. (2001'). We 
note that the data-set analysed by Pcrou x et al, (2Qll3) in- 
cludes that of lStorrie-Lombardi fc Wolfg ()2f)0fih . and a sim- 
ilar value for the D LA abundance was also reported by 
IStorrie-Lombardi fc Wolfe (2000). It is encouraging that the 
DLA abundances found in our simulations agree well with 
the observed range. 

As we discussed in Section r4.2l the DLA cross-section is 
falling off rapidly at Mtot ~ IQ^ h~^}AQ. Consequently, the 
DLA abundance per unit redshift can be read off from the 
cumulative abundance plot at Mtot = 10* /i~^Mq, provided 
that the underlying simulation can resolve this mass-scale 
well. For our highest resolution run at 2 = 3 (Q5), this is, in 
fact, the case. Here, the cumulative abundance has already 
flattened out at Aftot = lO*/i~^M0, so that a correction 
with equation ^ for a missed contribution by haloes on 
smaller mass-scales becomes unnecessary. 

In the left panel of Figure |S1 it is seen that the DLA 
abundance decreases as the wind strength increases from 
03 to P3, and to Q3, and as the resolution of the simula- 
tion increases from Q3 to Q4, and to Q5. This is due to 
the increasing slope of the power-law fits that we obtained 
in Section 0] As the slope of the power-law fit increases, 
the contribution from massive haloes becomes larger, while 
that from low-mass haloes becomes smaller. As a result, the 
cumulative abundance at the high-mass end of Figure |3 is 
largest for Q5, but when summed over all masses, Q5 ex- 
hibits the smallest total DLA abundance. 
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Figure 5. Cumulative abundance (or equivalently rate of incidence) of DLAs per unit redshift as a function of total halo mass. Left 
panel: Comparison of the results for the Q-series. The runs with a strong wind (Q3, Q4, Q5) have fewer DLAs than those with weak 
(P3) and no winds (OS). The transition from Q3 to Q4, and then to Q5, shows tha t high er numerical resolution tends to lead to fewer 
DLAs. The shaded region indicates the observed DLA abundance of lPeroux et al.l j22^). The solid slanted line shows the simulation 
result obtained bv lGardner et all 120011) . Right panel: Same as the left panel, but for the D- and the G-series, which all have strong 
winds. We also include the Q5 run as a thick solid line to ease comparison with the left panel. 



We also show the result from lGardner et alJ i200 j) as a 
solid slanted line in the left panel of Figure |K| They argued 
that their result would be consistent with the observations 
if the observed DLAs originate only from haloes with circu- 
lar velocities larger than Vc ~ 60kms~^ (which corresponds 
to Mtot « 2 X lO^'' h'^Mg). However, the good agreement 
between our improved simu lations and the o bserv ational de- 
terminations suggests that iGardner et alJ i200 J) probably 
overpredicted the DLA abundance due to the shallower slope 
they estimated for the relation between the halo mass and 
the DLA cross-section (see Section [4. H . 

In the right panel of Figure |S] the results for the D- 
and G-series are shown, with the values for Q5 included 
as a reference to ease comparison with the left panel. As 
the box-size increases from the Q- to the D-series, and then 
to the G-series, the resolution of the simulations severely 
degrades, causing an overprediction for the abundance, for 
the reasons discussed above. Therefore, we believe that Q5 
represents our most reliable estimate at z = 3 (leaving aside 
the question whether the feedback strength of the galactic 
wind adopted for the Q-runs is appropriate or not). 

5.2 Redshift evolution of DLA abundance 

In Figure [S] we show the evolution of DLA abundance from 
z = 4.5 to z = 0. In the left panel, the cumulative abun- 
dances of DLAs are shown as a function of total halo mass 
for redshifts z = 4.5, 3, 1, and 0. One can see that the con- 
tribution from massive haloes to the DLA population pro- 
gressively increases from high- to low-redshift as a result 
of the merging of haloes. In the right panel, we give the 
DLA abundance per unit redshift as a function of epoch. 



with values here read off at Mtot = 10^ h ^Mp, in the 
left panel. Observed data points from IPeroux et all ll200Jl 
and lRao fc Turnshek ( 2000) are also shown with error bars. 
Some of the exact simulation results are shown by the sym- 
bols, labelled with the names of the corresponding runs. The 
shaded region is our best-guess for a confidence region based 
on combining all of our simulation results. For reference, we 
show a power-law dN/dz = No{l + z)'' with A'o = 0.005 and 
7 = 2.5 as a long-dashed line, which describes the rate of 
evolution seen in the simulations well. It is encouraging that 
the above value of 7 is in good agreem ent with the observe d 
evolution of Lyman- limit systems of iPeroux et alJ i200ll) , 
where they report 7 = 2.45lg Q4. 

From z = 4.5 to 2 = 3, we see a decrease in the abun- 
dance by a factor of about two in both simulations and the 
observations, and the agreement between the two is very 
good, although the simulation points tend to fall slightly 
below the observations. From z = 3 to z — 1, the simulation 
(D5) suggests a further rapid decrease in DLA abundance 
by a factor of ~ 6, which is not seen at this level in the exist- 
ing observations. But the observational data at low redshift 
are still relatively uncertain, as indicated by the large error 
bars. The rapid decline is also reflected in the fact that f^Hi 
decreases from 0.66 {z = 3) to 0.14 (2 = 1) in D5 over this 
redshift range, a reduction of nearly a factor of 5 (see Fig- 
ure 0. If this significant decrease in the number of DLAs 
from 2 = 3 to 2 = 1 is real, it would partly explain why it 
is so difficult to find DLAs at 2 < 1. 

On the other hand, not much evolution is seen from 
2: = 1 to 2 = in the G5 simulation. This is related to the 
fact that JIhi in G5 does not decrease very much from 2 = 1 
(ilm = 0.19) to 2 = (Qm = 0.16). However, as discussed 
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Figure 6. Evolution of the DLA abundance from z = 4.5 to z = 0. Left panel: Cumulative DLA abundance as a function of total 
halo mass at redshifts z = 4.5,3,1 and 0. Rig ht panel: DLA abun dance per uni t redshift as a function of redshift. The data points 
with error bars are the observational data from lPeroux et alj i200l[) (crosses) and lRao fc TurnshelJ i200(]|) (open triangles a.t z < 1.5). 
The exact simulation results from some of the runs are indicated by the symbols with run-names. The shaded region is our best-guess 
for a confidence region based on combining all of our simulations. For reference, we show the power-law of dN/dz = No{l + z)^ with 
A^o = 0.005 and 7 = 2.5 as a long-dashed line. 



earlier, our power-law fits to the ctdla — A/tot relation are not 
well constrained for z = (and possibly for z = 1 as well), 
so the results ai z < 2 should be interpreted with caution. 
At 2 > 3, we saw that lower resolution runs tend to predict 
a larger abundance due to a shallower slope in the relation 
between the DLA cross-section and the halo mass, but it is 
not clear if other forms of systematic bias dominate at very 
low redshift for simulations with poor resolution. We will 
need yet higher resolution simulations with large box-sizes 
to make a more robust prediction of the DLA abundance at 
z < 2, and until then, it is not clear whether the current 
results for DLA abundance at z < 2, which tend to fall 
below the observational data, are trustworthy. This is why 
we have widened the shaded confidence region in Figure |S] 
significantly for z < 2. 



6 Hi column DENSITY DISTRIBUTION 
FUNCTION 

The column density distribution function f{N,X{z)) is 
defined such that f{N, X)dNdX is the number of ab- 
sorbers per sight line with Hi column densities in the in- 
terval [A^, A'' -I- dA'^] , and absorption distances in the interval 
[X, X -|-dX]. The absorption distance X{z) is given by 



x(.)^y^(i+.')^^d.'. (9) 

This definition is based on an argument by 

iBahcall fc Peebles! il969l) . who pointed out that the 
probability of absorption for a quasar sight-line in 
the redshift interval [z, z + dz] is dP oc (1 + z)^dr oc 



(1 -I- z)^[Ho/H{z)]dz = dX. In practice, if the comoving 
box-size of the simulation is AL, then the corresponding ab- 
sorption distance per sight-line is AX = {Ho/c){l + z)'^AL. 
For example, for AL — 10h~^ Mpc and 2 = 3, we have 
AX = 0.0534. 

Assuming that DLAs do not overlap along a sight-line 
through the simulation volume (which is a very good approx- 
imation given the small size of the simulation box, where 
the expected number of DLAs per sight-line at z = 3 for a 
10h~^ Mpc path is « 10~^), we can compute the Nm dis- 
tribution function by counting the number of grid-cells with 
column densities in the range [A'^, A'^ -I- dA']. In doing so, we 
are treating each grid-cell element as one line-of-sight. 



6.1 Hi column density distribution at 2 = 3 

In Figure |7] we show the Hi column density distribution 
function at 2 = 3. The solid triangles are the points directly 
measured from the s imulations. The ope n squares are the 
observational data of lPeroux et alJ i200ll for 2.7 < 2 < 3.5 
data), and the dashed curve is the fit to the same data based 
on a gamma-distribution: 



(10) 



The parameters of the fit are (/«, log A't , /3) = 
(0.0406, 21.18, 1.10) iPeroux et aP J200ll for 2.7 < 2 < 3.5 
data). We note that all data by IStorrie-Lombardi fc Wolfd 
(2000) are included in that of Peroux et al.'s. 

In the panel for 'Q3' (upper right corner), we also show 
the result of different smoothing methods, using crosses 
(uniform cloud-in-cell distribution with I = [47r/3]^''^s) 
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Figure 7. Hi column density distribution func tion at 2 = 3. The so lid triangles are the points measured directly from the simulations. 
The open squares are the observational data of lPeroux et al ] <200ll for 2.7 < 2 < 3.5 data), and the dashed line is the fit to the same 
data based on a gamma-distribution. In the panel for 'Q3', the results with different smoothing methods are shown in crosses and open 
triangles. See text for details. 



and open triangles (uniform clouds-in-cell distribution with 
£ — |[47r/3]^'''^s). The former method (crosses) results in 
higher values of f{N) at lower column densities because 
it smoothes the gas mass into broader regions. The SPH 
smoothing method agrees with the latter calculation method 
(open triangles) better. 

The agreement between the observations and the simu- 
lations Q3, Q4, Q5, & D5 at logA''Hi > 21 is generally very 
good. Results from runs of increasing resolution (Q3, Q4, 
and Q5) are consistent with each other to a high degree. 
The run with no wind (OS) somewhat overpredicts the dis- 
tribution function at large A^hi values, but as the galactic 
wind strength increases from 03 to P3, and then to Q3, 
the high column density systems become less abundant and 
the agreement between the simulation and observations im- 
proves. At intermediate column densities (20 < logA^Hi < 
21), it seems that the simulated distribution function falls 
short of the observational estimate. Given the consistent be- 
haviour in Q3, Q4, and Q5, our result appears not to be 
affected by resolution, although this cannot be completely 
excluded. We will discuss this point further in the next sub- 
section, when we consider the data at z = 4.5. It is clear 



however that G4 and G5 do not have sufficient resolution at 
z = 3 to resolve DLAs. 

6.2 Hi column density distribution at z = 4.5 

In Figure |H1 we show the Hi column density distribution 
function at 2 = 4.5. As before, the solid triangles are the 
points measured in the si mulations, and the op en squares are 
the observational data of lPeroux et al.l ll200ll . for 3.5 < z < 
4.99 data). The long-dashed line is the gamma fit to the same 
observational data of 3.5 < z < 4.5, and the short-dashed 
line is the fit to the data for 2.7 < z < 3.5 for reference. 
The values of the fit parameters for the 3.5 < z < 4.5 data 

is {f:,,\ogN,,l3) = (0.2506, 20.46,0.80). 

Observational studies lIStorrie-Lombardi fc Wol^l2000l : 
iPeroux et al]l200ll) indicate that there are fewer high A'hi 
systems (log A'hi > 21) at z > 3.5 compared with 2.7 < 
z < 3.5, and that the distribution function becomes steeper 
at 2 > 3.5. However, we do not see such a reduction of 
high A^Hi systems in our simulations from z = 3 to z = 
4. In fact, the highest resolution simulation in our series 
(Q5) suggests that f{N) is slightly higher (but steeper at 
the same time) at z — 4.5 compared to 2 = 3. Note that 
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Figure 8. Hi column density distribution function at z = 4.5. The solid triangles are the points measured in the simulations, and the 
open squares are the observational data of lPeroux et al ] <200ll for 3.5 < ^ < 4.99 data). The long-dashed line is the gamma fit to the 
same observational data of 3.5 < z < 4.5. For comparison, the short-dashed line is the fit to the 2.7 < z < 3.5 data. 



the agreement between runs with different resolution (Q3, 
Q4, Q5, R3, and R4) is impressive, showing that the results 
are well-converged. The degree of the increase in f{N) from 
2: = 3toz = 4.5is somewhat larger for the intermediate A'hi 
systems (20 < logA^ni < 21), leading to a slight steepening 
of the overall f{N). One interpretation is that observational 
studies have not yet discovered sufficient numbers of DLAs 
with very high column density (log A^hi > 21) to accurately 
estimate the evolution of f{N) at z > 3, because such high 
A^Hi systems are intrinsically rare and therefore difficult to 
find. 

The evolution from z = 4.5 to z = 3 is presumably 
driven by the combined effect of the depletion of the abun- 
dance of small haloes by merging, the accumulation of feed- 
back by galactic winds, and the reduction of cooling effi- 
ciency. Small haloes of comparatively low column densities 
(logA''Hi < 20.5) merge into larger systems from z = 4.5 to 
2 = 3, forming more massive systems with higher column 
densities, as the hierarchical structure formation scenario 
suggests. However, strong feedback by galactic winds ejects 
neutral gas from the star-forming regions, thereby acting to 
reduce the column densities of all systems. In addition, the 
efhciency of cooling rapidly decreases towards lower redshift, 



reducing the rate at which gas can cool out of haloes and 
become neutral. 



6.3 Hi column density distribution at lower 
redshift 

In Figure |^ we show the Hi column density distribution 
function at z = 1 and z — 0. Again, solid triangles 
give the points measured in the simulations, and the open 
squares are the observational data of Pcroux ct al. (20o3, 
for 0.008 < z < 2.0 data). The dashed curve is the gamma- 
distribution fit to the same observational data with param- 
eters (/.,log A.,/3) = (0.0870,20.76,0.74). For comparison, 
we also show the observational gamma fit for 2.7 < z < 3.5 
data as dotted lines. 

As we pointed out earlier, the results of the G-series 
are substantially affected by resolution, causing an under- 
estimate of the DLA abundance if it is not corrected us- 
ing equation (If the abundance is corrected with an 
ill-determined relation between the DLA cross-section and 
the halo mass, then low resolution can also lead to an 
overestimate of the abundance, as we discussed in Sec- 
tion|3) With the mass resolution of the G-series, haloes with 
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Figure 9. Hi column density distribution function at 2 = 1 
and 2 = 0. The solid triangles are the points measured in the 
simulations, and th e open squares are the observational data of 
IPeroux et all i200ll. for 0.008 < z < 2.0 data). The dashed curve 
is the gamma-distribution fit to the same observational data. 
For comparison, we also show the observational gamma fit for 
2.7 < z < 3.5 data as dotted lines. 

Mtot < 10"' h-^M@ are not resolved, and a significant frac- 
tion of the DLA population is therefore missed. Curiously, it 
is seen that the column density distribution function in fact 
rises from z = 1 to z = in the G-series. This is probably 
because haloes with masses above the mass resolution limit 
are only beginning to form at these low redshifts as the non- 
linear mass-scale increases with decreasing redshift, and the 
simulation is finally "catching up" to reproduce the neutral 
gas content in these haloes. 



7 DISCUSSION 

We have used state-of-the-art hydrodynamic simulations of 
structure formation to investigate the abundance of DLAs 
in a ACDM universe. Our study represents a first attempt to 
apply a large series of simulations to this problem, probing 
an unprecedented range in both mass and spatial scales, 
enabling us to quantify systematic effects due to numer- 
ical resolution. Furthermore, we improved the simulation 
methodology by adopting a novel formulation of SPH (see 
ISprineel fc Hernauistl20oil that minimises systematic inac- 
curacies in simulations with cooling, and by using an im- 
proved model for the treatment of the multiphase structure 



of the ISM in the context o f star formation and feedback 
JSprineel fc Hernauistll2003all . 

By comparing our results for DLA abundance in a series 
of runs as a function of resolution and feedback strength, we 
were able to demonstrate that insufficient resolution, or a 
lack of a proper treatment of effective feedback processes, 
leads to an incorrect estimate of the relation between the 
DLA cross-section and the halo mass. This likely led to an 
overestimate of the DLA abundance in earlier studies, for 
the reasons we discussed in det ail in Section f4.1l 

IProchaska fc Wolfd J200ll) pointed out that the ob- 
served velocity width distribution cannot be reproduced if 
the relation between the DLA cross-section and the halo 
mass derived from these earlier SPH simulations is used. 
They also suggested that one possibility to remedy this in- 
consistency was to suppose that the relationship between the 
DLA cross-section and the halo mass was incorrectly deter- 
mined. This is exactly what we find in our current study. 
If we use the new relation found in our highest resolution 
simulation, we obtain a DLA abundance that is consistent 
with observations, which is very encouraging. The slope of 
the relation that we infer from our simulations at 2 = 3 is 
in t he range of 0.7 — 1.0, which co i ncides with the range 
that iHaehnelt. Steinmetz. fc RauchI l|2QQCfl derived by re- 
quiring their model prediction of velocity width distribution 
to match the observed one. 

However, while our simulations reproduce the DLA 
abundance a.t z — 3 very well, our predictions at 2: < 2 
are not equally reliable because they are based on simula- 
tions with larger box-sizes and lower resolution. This is also 
evident from the poor agreement between the simulated and 
the observed Hi column density distribution function at low 
redshift. To make the predictions at low-redshift more ro- 
bust will require large box-size simulations with yet more 
particles. 

The fact that we were able to reproduce the observed 
DLA abundance very well at redshift z — 3 has significant 
implications for the cold dark matter model. This result sug- 
gests that DLAs arise naturally in a ACDM universe from 
radiatively cooled gas in dark matter haloes with correct 
abundance. This is related to the "substructure problem" 
that is posed against the CDM model, based on the no- 
tion that the number of satellite galaxies observed around 
the Milky Way seems to be much smaller than, and hence 
in contradiction to, the large number of dark matter sub- 
structures predicte d by high-resolution N-body simulations 
jMoore et alJIlQQSt). 

IstoehTwhite. Tormen fc Springell ^2002^ have shown 
that the seriousness of this discrepancy was overstated ini- 
tially, but the high abundance of low mass haloes and dark 
satellites in CDM models remains puzzling, given for ex- 
ample the shallow faint-end slope of the galaxy luminos- 
ity function. Therefore, many models have been invoked 
to suppress the formation of dwarf satellite galaxies, in- 
cludi ng supernova feedback that ejects gas jDekel fc Silkl 
Il98(fl . reheating of the intergalactic medium w hich sup- 
presses gas infall jBullock. Kravtsov. fc Weinb erg 2001J), 
and photo ionisation of gas by UV b ackground radiation 
jThoul fc Weinberg 199 6: Q uinn. Katz. fc Efstathiou 199^. 
Full cosmological hydro-simulations have also been used to 
argue that feedback effects by UV background radiation 
and supernovae may account for the suppression of the for- 
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mation of low-mass galaxies rel ative to the steeply risinR 
dark matter halo mass functio n iChiu. Ostriker. fc Gnedir] 
I2OOII : iNaeamine et al.l 1200 llj) . While we cannot make a 
strong statement on this problem at low-redshift, our re- 
sults suggest that the CDM model does not have dif- 
ficulty at z = 3 with respect to the number of neu- 
tral gas clumps from which stars form. Considering 
that there is in general good agreement between ob- 
servations and theore tical studies of Lyman - break galax- 
ies at z = 3 (e.g. IMo Fukugital Il996t iBaugh et alJ 
19981: IJing fc Sutdll998l:lKatz. Hernauist fc Weinberdll999t 
Kauffmann et al." 1999^ 'Mo. Mao, fc Whitel ll999l: iNagaming 
2002; Weinberg, Hernauist & Katz 2002), the ACDM model 
seems to be in a very good shape at z = 3. 

Another interesting result of our study is that the 
break in the relation between DLA cross-section and halo 
mass seems to occur at a mass-scale of Mtot ~ 10** — 
lO*'^/i~^M0. DLAs do not exist in haloes below this 
mass-scale at z = 3 and z — 4.5. Because we mea- 
sured this effect consistently in both the 'Q5' and 'R4'- 
runs, which have sufBcient resolution to describe haloes 
with Mtot ~ 10* h~^MQ well, it is clear that this is not 
caused by a resolution effect. Note that this mass-scale, 
which corresponds to circular velocities of 10 — 15kms~^ 
and virial temperatures of about 10* K or even slightly 
lower, is smaller than sug g ested by th e earlier studies of 
iQuinn. Katz. fc Efstathioul Jl996j) and [Th oul fc Weinberg 
Jl99d) rSee also IWeinberg. Hernauist fc Kata H997l) l. The 
physical origin of the break could lie in both the radiative 
processes of cooling and UV-heating, and in supernovae feed- 
back, but we expect that the former effects dominate. This 
is because feedback will proceed only if at least some neu- 
tral gas is built up, so that star formation can occur at some 
(low) level. Only after that can feedback quench the further 
accumulation of neutral gas, and possibly also affect neigh- 
bouring haloes. Therefore, while one can easily understand 
why supernovae feedback reduces the amount of neutral gas 
in small halos, it is not clear how it could lead to the com- 
plete absence of neutral gas below a certain mass-scale. On 
the other hand, a relatively sharp break can be expected 
from the properties of the cooling function, and the heat- 
ing of the gas due to the UV background. Both processes 
depend very sensitively on the virial temperature of halos 
in the range 10* — 10^ K. Further studies will be needed to 
disentangle the relative importance of the various physical 
effects in this regime more precisely. 

Assuming that feedback effects due to star formation 
nevertheless influence the mass-scale of the break, we might 
have underestimated it if the feedback we employed in our 
simulations was not strong enough. However, we saw in 
our results of Figure that the current simulations with 
the 'strong wind' model already seem to underpredict Qui 
slightly, therefore even stronger feedback would make the 
neutral gas fraction in the Q-series even smaller and exac- 
erbate the discrepancy with the observational estimate of 
Qm at z — 3. Also, there is only very limited room to in- 
crease the strength of the UV background flux, because for 
our chosen normalisation the mean opacity of the Lyman-a 
forest is in good agreement with observations. Another way 
of constraining the feedback strength is to compare the sim- 
ulated galaxy luminosity function with the observed one, to 
see if the flatness of the faint-end of the luminosity function 



can be reproduced. A preliminary analysis suggests that our 
present simulations do not yet yield a luminosity function 
with a faint-end as flat as observed at low-redshift; therefore 
the feedback model in our current simulation might not be 
adequate to correctly account for the galaxy population at 
z = 0. The results of this analysis will be presented else- 
where. A solution to this problem will help to further con- 
strain the physical processes that suppress the DLA cross- 
section in low-mass haloes. 

In the present paper, we focused on the abundance of 
DLAs, without carrying out a detailed study of their inter- 
nal structure. Therefore, we are not able to completely rule 
out the possibility that the DLAs we find in our simulations 
are in fact originating from very small disks at the center 
of dark matter haloes. Our present work is complementary 
to that of Haehnelt et al. (1998), in the sense that we are 
able to sample a large statistical ensemble of absorbers in 
comoving volume of > {10h~^ Mpc)'', while they studied a 
smaller number of objects in greater detail. Nevertheless, 
our highest resolution run (Q5) with 10h~^ Mpc box has 
a spatial resolution of 1.2h^^ kpc which is comparable to 
that of Haehnelt et al. (1998), and the results from vary- 
ing resolution runs (Q3, Q4, Q5) agree with each other well. 
Therefore, given that our simulations are based on the CDM 
model, our results support the idea that DLAs arise from ra- 
diatively cooled protogalactic gas clumps embedded in dark 
matter haloes. While our spatial resolution is approaching 
comoving sub-kpc scales at z = 3, it is probably not yet suf- 
ficient to accurately resolve the disk structures embedded in 
the central few kpc regions of dark matter haloes. Kinemat- 
ical studies of the internal structure of DLAs will require yet 
higher resolution simulations than those used in the present 
paper, and possibly a more sophisticated physical model for 
the ISM and star formation. 
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